A reaction-diffusion model for the hydration of calciumsulphate (gypsum) and 

microstructure percolation. 



F. Tzschichholz 1 and H. J. Herrmann 1 

1 Institut fur Computeranwendungen I, Universitat Stuttgart, 
Pfaffenwaldring 27, D-70569 Stuttgart, Germany 
(February 1, 2008) 

We have numerically investigated a reaction-diffusion model for the hydration of calcium- 
sulphate (gypsum). The simulations were conducted for two and three dimensional systems. 
While the dissolution of anhydrous gypsum is considered irreversible at a finite rate the pre- 
cipitation/dissolution reaction for the calciumdihydrate is considered reversible. The latter 
reaction is assumed to be controlled by the dihydrate's equilibrium solubility and the abillity 
of the system to react on supersaturation only at a certain velocity described by the reaction 
rate constant of precipitation. For d = 2 we find at early times an accelerated hydration period 
followed by a maximum and a decreasing hydration rate. For large times the ionic product of 
involved species assumes closely the value of the di-hydrate equilibrium solubility. Calculated 
model micro-structures exhibit typical features such as inner and outer hydrate products, in- 
duction and dormant period as well as bridging. Furthermore we find that the overall chemical 
reactivity as a function of initial anhydrous (volume) concentration p exhibits a maximum 
close to the percolation point of the underlying lattice. Employing a rescaling procedure we 
find two percolation thresholds in d = 2, p™ in = 0.44 ± 0.015 and p™ ax = 0.77 ± 0.02, for the 
initial anhydrous gypsum concentration between which percolating dihydrate structures can be 
attained. For d = 3 we find p™ in = 0.10 ± 0.02 and p™ ax = 0.95 ± 0.02. 

PACS number(s): 61.43.-j, 81.35.+k, 81.30.Mh 



I. INTRODUCTION 

Gypsum is perhaps one of the oldest crafting and 
building materials human kind has cultivated beside 
wood, iron, and stone. There exists a long tradition 
of using this material in arts, medicine, paleontology, 
archeology, as well as a light-weighted building ma- 
terial. Its widespread use is likely due to its natural 
abundance, its almost total flexibility in applications, 
its chemical inertness and certainly its low costs. For 
an overview of the various applications we refer the 
interested reader to the survey of Wirsching [1] 
Gypsum belongs to the group of calcium-sulphates 
with its various hydrates. A material class exhibiting 
some common phenomenological properties are the 
calcium-silicates, mostly termed as cementious mate- 
rials. Cementious materials are, however, of much 
higher intrinsic strength than gypsum due to their 
chemical shrinkage and are therefore in technical ap- 
plications mostly prefered to gypsum. While most 
cementious materials do precipitate as a gelantineous 
phase, calcium-sulphates appear always in crystalline 
modifications (needles). 

Both calcium-sulphates as well as calcium-silicates are 
already well characterized on a molecular level [2,3]. 
However, it appears to us that gypsum is less un- 
derstood compared to calcium-silicates on the macro- 
scopic scale. For more recent experimental investi- 
gations we refer the interested reader to Ref. 4. For 
cementious materials microstructure computer mod- 
els have been recently devellopped [5-8]. Though 
the cementious systems are much more complex from 
a chemical point of view their huge economical im- 
portance as the building material of choice has fo- 
cussed the attention on the modeller's side. In the 



on diffusion-reaction calculations. In the practice, the 
most important feature of the setting of gypsum is 
the formation of rigid structure of the hydrates. The 
necessary condition for this to occur is a connected 
structure, in other words, that the hydrate aggregate 
percolates. In this paper we will therefore investigate 
under which conditions (concentration of a hydrate) 
the end product of the reaction forms a percolating 
cluster. 



II. PHYSICOCHEMICAL ASPECTS 

The hydration of anhydrous calciumsulphate, CaSO^ 
- in the following CS, is based on dissolu- 
tion/precipitation reactions forming various calcium- 
sulfohydrates [1]. The hydrates, CaSO^ ■ nH20, are 
distinguished by the molar amount n of bound crystal- 
water. In connection with the amount of physically 
bound water one observes different physico-chemical 
properties as for example crystal-symmetries, densi- 
ties, or solubilities. The two most important hydrates 
are the semihydrate, n = 1/2, and the dihydrate, 
n = 2. There are two 'modifications' of the semi- 
hydrate called the a— and (3— semihydrate. On the 
side of the anhydrous calciumsulphates three forms 
are distinguished: soluable CS (A III), dead CS (A 
II), and high temperature CS (A I). Typical solubil- 
ities at room temperature and normal pressure are 
2 gramm / 'liter for the dihydrate, 2.5 gramm/ liter 
for All, 5.8 gramm /liter for a-semihydrate, and 
7 .6 gramm / liter for /3-semihydrate. The solubili- 
ties for a— ,f3— semihydrate and All decrease mono- 
tonically with increasing temperature while the di- 
hydrate exhibits a flat maximum at temperatures 
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which consists almost only of dihydrate, b) anhydrous 
stone consisting almost only of A II, and c) technolog- 
ical mixing forms consisting mostly of A II and A III 
and semihydrates. 

Stochiomctrically one has the following net reac- 
tions: 

hydration of anhydride to dihydrate 



CaSOi + 2H 2 0^ CaS0 4 ■ 2H 2 0, 
hydration of semihydrate to dihydrate 



(2.1) 



2CaS0 4 ■ \ H 20 + 3H 2 -» 2CaSO A ■ 2H 2 0. (2.2) 

The implementation of the foregoing reactions in a 
reaction-diffusion model requires explicit knowledge of 
the dissolution/precipitation reactions, i.e., their ki- 
netic rate equations. Observing that these rate equa- 
tions are not explicitely known to us we make the 
following assumptions. 

1. Solubilities. The present calciumsulphate sys- 
tem is strongly electrochemical. In the following 
we will assume nevertheless that all electrostatic 
interactions between solvated ions can be dis- 
regarded. Solubilities are therefore given by 
ionic products of concentrations. In particular, 

S cZso 4 ~ 10 3 Gramm/ liter > S Ca sOi-2 h 2 o ~ 
2Gramm /liter, S CaSOi .i H20 « Q,hGramm /liter. 
All values are equilibrium values at room tem- 
perature and normal pressure. The value for 
ScaS0 4 ^ S supposed to correspond to an in- 
finitely high solubility. 

2. Diffusion constants. Correspondingly we em- 
ploy diffusion constants as observed for 'infinite 
dilution' neglecting all possible electrochemi- 
cal influences. This is a simplification for liq- 



uid electrolytes. In particular, D L 
1CT 10 m 2 /s, D so l~ « 10.7- 1(T 10 m 2 /s. 



7.9 



Solvation. Considering the process of solvation 
it becomes clear that the above hydration reac- 
tion needs to be split up at least into a disso- 
lution and a precipitation reaction. The formed 
ions are solvated ions binding a certain amount 
of water. Just let us write down the dissolution 
reaction, 



The numbers vc a and vso 4 do characterize 
the formed hydrate shells of the calcium- and 
sulphate-Ions. One has vca — 6. We are cur- 
rently not aware of the value for vso A - For the 
simulations we tentatively assume vsOi — 6. 

4. Specific molecular volumes. The specific 
volumina of the involved solid and liquid 
phases under normal conditions were ob- 
tained from the literature [1], in particular: 
52.7 • 10 liter /mol, vc a so 4 - 



VCaSOi 

74.1 • 10- 3 liter/mol, 



v. 



CaSOi-0.5H 2 O 



10 3 liter/mol, v, 
10~ 3 liter/mol and vr 2 o = 18 • 10" 



(/») 

CaSO 4 0.5H 2 O 



■2H 2 = 

55,7 
= 52.7 
3 liter /mol. 



5. Reaction constants. The reaction constants were 
chosen such that firstly the back reaction of the 
anhydrous dissolution becomes negligibly small 
(irreversible reaction) and secondly the precipi- 
tation/dissolution of dihydrate is adjusted to the 
experimentally observed value for dihydrate sol- 
ubility. For calciumsilicates typical scales for the 
involved surface reaction rates are 10~ 6 m , . 

raoL-s 

We have tried to orientate ourselves on these 
magnitudes. 



III. THE MODEL 

In the following we consider only the hydration of 
anhydrous gypsum towards dihydrate. As precise re- 
action mechanisms and rates are not known to us it 
appears reasonable to model just one dissolution and 
one precipitation reaction, i.e., Eqs. (2.3) and (2.5). 
The general model setup has been described elsewhere 
[8] . For the presented calculations we employed a time 
integration step At = 10 _1 s and a spatial resolution 
Ax = 10~ 4 to. The calculations were performed in 
d = 2 for a system of size 100 x 100 and in d = 3 
for sizes of 50 x 50 x 50. The dissolution rate con- 
stant was prescribed and fixed at kdi SS = 10~ 3 
(absolute scale) for a somewhat lower precipitation- 
rate-constant kp rec = 10~ 2 • kdiss- The by a factor 
100 lower precipitation rate constant was choosen in 
order to mimick the experimental observation, of a 
much faster dissolution than precipitation. 



{CaS0 4 ) {s) + (v C a + vsojH 2 {e) -> 
(Ca 2+ , u C aH 2 0) (aq) + (SOl~ ,v SOi H 2 0) {aq) , 

und the precipitation reaction, 

(Ca 2+ ,v Ca H 2 0) (aq) + {S0 2 -,usoH 2 0) (aq) - 
(CaS0 4 ■ 2H 2 0) (S) + [y Ca + v SOi - 2)H 2 (l) . 



(2.3) 
(2.4) 



IV. RESULTS 



A. Results for d = 2 



Fig. 1 represents the volume concentrations for an- 
(2-^ydrous gypsum, dihydrate, and ions (left to right). 
(2.6fJ,ed means high - blue low concentrations. 
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The cases k<p T&c — 0.1 k(H SS Hind k prec — k d -i ss 
hibit similarities where the early hydration undergoes 
an 'induction' period. In this induction period the 
ion concentrations rise steeply due to a relatively high 
anhydrous gypsum dissolution. After passing a max- 
imum ion concentration the main reaction period be- 
gins to emerge in which 'outer' dihydrate starts to 
precipitate blobking further reactions of the anhydrid. 
Correspondingly ion concentrations drop to an almost 
stationary value, considerably slowing down the whole 
further hydration in the vicinity of the anhydrous gyp- 
sum grains. Such behavior is also known from cal- 
ciumsilicatcs during the formation of inner hydrates 
(experimentally and numerically). The for k prec over 
two decades almost constant ion concentration cor- 
responds essentially to the equilibrium solubility for 



drous gypsum is precipitated relatively quickly again. 
The velocity with which the precipitation reaction re- 
acts with respect to the supersaturation is controlled 
by the rate constant k prec . 



In order to study the influence of the initial volume 
of anhydrous gypsum p on the hydration process a set 
of calculations was conducted. The calculations were 
performed for k dlss = lCT 3 ^ and k prec = 10^ 2 k diss . 
The basic investigated question was for which p criti- 
cal percolation of the hvdrated vhase in the long run 
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FIG. 4. Volume densities after 8 h of simulated hydra- 
tion. k di3S = 1(T 3 ^, k prec = W~ 2 k dias , L = 100, anhy- 
drate (o), water (A), dihydrate (□). All other parameters 
as in Fig. 1. 

As a first investigation the total volume densities 
after eight hours of hydration time were measured in 
a range between p = 0.4 and p = 0.8, see Figure 4. 
The system size was 100 x 100. The numerical effort 
per data point was about 40 minutes cpu time on a 
2GHz PC. 

For initial anhydrous gypsum densities less than 
p = 0.54 the final remaining gypsum is less than 
0.02, i.e., the dissolution reaction is complete. Above 
p = 0.6 one observes a strongly increasing fraction 
of non dissolved gypsum, i.e., the reactions become 
increasingly incomplete. As water is not treated as 



limiting rcactant per se in the calculations there must 
exist a non-stochiometric reason for the incomplete 
reaction. The amount of precipitated dihydrate (□) 
confirms this. The corresponding curve undergoes a 
maximum in reactivity for p values between 0.6 and 
0.65 although sufficient amounts of water are system 
wide available. This effect is most likely due to a 
screening of the initial anhydrous gypsum. 
Above p c = 0.592 increasingly parts of the chemi- 
cally reactive interface between water and gypsum are 
blocked by the gypsum itself. For p > p c one has to ex- 
pect a decreasing specific surface area for the gypsum 
- which corresponds directly to a decreasing chemical 
reactivity. Let p denote the total volume fraction of 
initial anhydrous gypsum. As the hydrate density is 
a more or less continous quantity p(x,p) € [0, 1] one 
has to term the statement of percolation somewhat 
different than 'usual'. Looking at the densities p, we 
define a function R(x,p,r) with the property, 



R(x,p,r) = @(p(x,p) - r) 



(4.1) 



(with 9 being the Heaviside step function). The 
function R(x,p,r) maps just all those volume ele- 
ments to unity (red) for which the corresponding 
hydrate density exceeds a certain amount r, other- 
wise it becomes zero (blue). The following images 
show the hydrate pixels R(x, r) at r = 0.7 for p = 0.4, 
p = 0.5, p = 0.7, and p — 0.8 (from left to right). 




FIG. 5. Hydrate pixels at r — 0.7 for initial anhydrate volume concentrations p — u.4, p = u.o, p = V. I, and p 
(from left to right). A red pixel corresponds to an asymptotic dihydrate volume concentration larger than r whereas 
blue pixel represents a concentration lower than r. Parameters as in Fig. 4. 



From these images it is apparent that for very high 
gypsum densities, i.e., p = 0.8, the hydrate does 
not percolate that easily as compared for example to 
p = 0.6. This indicates that at least two percolation 
transitions for the hydrated phase do exist. 

We employ volume concentrations in this work. Ob- 
viously the solid phase concentrations define a mea- 
sure for the distinct properties. How to define crite- 
ria in terms of phase concentration(s) that allow an 
analogy to percolation phenomena? In conventional 
(site-)percolation each volume element is either com- 
pletely occupied by species A or B at random in an 
uncorrelated way with fixed occupation probabilities 
v a and 1 — D4 respectively. The constant occupation 



tion exhibits scale invariance (sclfsimilarity) [9]. Vice 
versa one could employ this scale invariance as a tool 
in order to estimate the critical percolation probabil- 
ity pa- We will follow this approach in the present 
paper. Because the hydrate concentrations directly 
correspond to non-constant occupation probabilities 
one could expect correlated percolation. 

The microstructural and topological properties of 
the transformed phase strongly influence overall (vol- 
ume) properties such as elastic properties. If -as in our 
case- the unreacted phase owns a small or even van- 
ishing elastic modulus (anhydrous gypsum is usually 
a powder) then the percolating phase controls entirely 
overall properties. Thereforethe percolation transition 
of the spatial dihydrate distribution is of central im- 
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question arises on how to estimate and characterize a 
possible percolation transition as for our case. We in- 
troduce the typical volume density r above which all 
hydratcd volume elements are regarded as occupied 
sites in the sense of percolation. Let us denote k(p, r) 
as the true overall density of such sites. In general 
k(p, r) and r will have different values because the 
overall density k depends on the anhydrous volume 
fraction p while r is a variable. However, if we guar- 
antee by some algorithm that the occupied sites do 
always percolate for arbitrary values of r (scales) by 
adjusting the initial p, then the overall density k must 
equal the typical density r. In other words k = r 
which is a fixed point of the map k{p, r). 

In order to investigate this question quantita- 
tively more extensive calculations have been per- 
formed. One can calculate from the pixel-distribution 
R(x,p, r) two main quantities of interest a) its overall 
density 



fc(P,r) = \ J d 3 xR(x,p,r) (4.2) 



and b) the information whether the pixel- 
distribution R is percolating itself (via a burning al- 
gorithm [9] ) . The set of rules without the fixed point 
condition is the following: 




Anhydride fraction p 
FIG. 6. Hydrate pixel volume density,fc(p, r), versus ini- 
tial anhydrate concentrations, p, for three different cutoffs 
r = 0.5(o), r = 0.593(+) and r = 0.7(A). Parameters as 
in Fig. 4. 



The data points link the dependence of the density 
k and initial anhydrous volume density p for various 
cutoffs r under the condition that R is percolating. 

It is more meaningful to average all points for the 
same r belonging to a percolating cluster, and calcu- 
late their mean and variance values, which results in 
one single point with error bars per value of r. This 
is shown in Fig. 7. 



1. Choose an initial 



0.593, 



1 



2. Choose a sufficient low initial gypsum density, 
i.e., p = 0.3 in order to find later On Pmin- 



3. Do the simulation until stationary concentration 
fields are reached ( here t — 3 • I0 4 s). 



4. Calculate R(x,p 7 r) and k(p,r). Determine 
whether R(x,p,r) is percolating or not, and 
choose accordingly a decremented/incremented 
p. A hyperbolic decreasing interval incre- 
ment/decrement was chosen. 



5. Continue with 3. until the increment interval 
becomes lower than 5p = 10~ 3 . 



^0.8 
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FIG. 7. Averaged hydrate pixel volume density, k(p c , r), 
versus cutoffs r for the various critical anhydride densities 
as observed from Fig. 6. Points that intersect the line 
k(p c ,r) = r define a critical value k c (p c ,k c ). Parameters 
as in Fig. 4. 



6. Choose a different initial value for r, and con- 
tinue with 2. 

The algorithm works similar for p max . 

Fig. 6 shows some tvmcal plots of hvdrate densities 



We now impose the condition k(r 7 p c ) = r on the 
data points shown in Fig. 7. We find a common 
critical value k c = 0.62 (+0.03/ - 0.04). From the 
lowest and highest value for k c one obtains the de- 
sired values for the minimum and maximum perco- 
lation threshold p™" = 0.44 (+0.015/ - 0.015) and 
raax = qj 7 ( +om / _ 0.02). This means in narticu- 
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In the foregoing paragraph we investigated a few 
percolation properties in d = 2 depending on the ini- 
tial anhydrous gypsum concentration. We found that 
for the case of a diluted (p < p™ m ) and dense ini- 
tial anhydride packing (p > p™ ax ) no percolation of 
the hydrated phase will occur, i.e. the hydrated gyp- 
sum becomes useless for most technical applications. 
We investigated this problem also in d = 3. To this 
end the reaction-diffusion model and the burning al- 
gorithm were modified. 
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Anhydride fraction p 
FIG. 9. Hydrate pixel volume density, A; (p, r), versus ini- 
tial anhydride concentrations, p, for seven different cutoffs 
r = 0.25(*), r = 0.27(o), r = 0.29(V), r = 0.312(+), 
r = 0.33(A), r = 0.35(D) and r = 0.37(0). Other pa- 
rameters: k diss = 10 _3 ^||, kprec = 10~' 2 k dtss , L — 50 in 
d = 3. 
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FIG. 10. Typical volume density, k(p c ,r), plotted ver- 
sus r. Points that intersect the dashed line k(p c , r) = r 
define the critical k c — 0.25 ± 0.02 in d = 3. Parameters 
as in Fig. 9. 



FIG. 8. A typical dihydrate backbone in d — 3, system- 
size 20 3 . 



In Fig. 8 we show a small hydrate backbone, i.e. 
sites not belonging to a spanning cluster have been 
removed. One clearly sees the holes in the microstruc- 
ture constituting the pore space. In order to deter- 
mine the further percolation properties we employed 
the same numerical procedure and parameters as ex- 
plained in Sec. IV A with the exception of the system 
size being now 50 3 . 

The numerical effort was about 7 davs CPU time 



The scatter of data for each r value appears much 
smaller than in Fig. 6. Of course are the limits for 
k(r,p — > p c ) smaller now, see Fig. 6, because the 
thresholds for uncorrelated percolation in d — 3 are 
much smaller than in d = 2. In d = 3 there exist 
much more connected pathways for ion diffusion and 
solvation. Therefore it is reasonable to expect that a 
possible upper percolation threshold p" lax is shifted 
to higher anhydrous concentrations as compared to 
Fig. 6. Again in Fig. 10 we show mean and variance 
values for k versus r. 

Imposing the condition k{r,p c ) = r we obtain 
within the error bars k c — 0.25 ± 0.02 implying 
p min = 010 ± o.02 and p™ ax = 0.95 ± 0.02. As p™ ax 
is verv high compared to d = 2 this could nlav a cer- 
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V. CONCLUSION 

We studied percolation properties of a reaction- 
diffusion model for the hydration of calciumsulphate. 
From a kinetic point of view the basic effects of the dis- 
solving and precipitating microstructures are a) that 
they control the ion diffusion transport to a large ex- 
tent and b) that they form a spanning backbone from 
one systems end to the other, thus allowing the tyrans- 
port of momentum/stress. 

We expect that for p < p™ ln and p > p™ ax mi- 
crostructures are unable to carry stress at all. It 
should be noted that for d = 2 and d = 3 p™ ax is 
not reachable employing mono-disperse packings. 

Our results also show that there exists an optimum 
initial anhydrous concentration for which the total 
amount of dihydrate reaches a maximum, i.e. shows 
maximum overall chemical reactivity, and possibly the 
best mechanical properties. 

The relatively high computing time for d = 3 indi- 
cates that one possibly should consider in the future 
a different evaluation method than exact enumeration 
in order to solve the reaction-diffusion equations. 
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